// pt20m.C -- post ta pid going on -- z = 20, merging rootfiles using TChain

// p[0]: aoz; x[0]: beta
// return p/Z = poz(beta, aoz);
double pv(double *x, double *p);
double pbg(double *x, double *p);

void pt20m(){

  #include "/home/asia/Lab/etfana/etfana/config/apr2023/cmd/cut/ca21.C"
  #include "/home/asia/Lab/etfana/etfana/config/apr2023/cmd/cut/ca22.C"
  #include "/home/asia/Lab/etfana/etfana/config/apr2023/cmd/cut/ca23.C"
  #include "/home/asia/Lab/etfana/etfana/config/apr2023/cmd/cut/ca24.C"

  // TFile *f = new TFile(getenv("TMP1"), "UPDATE");
  // TTree *etf = (TTree *)f->Get("etf");
  TChain *etf = new TChain("etf");
  etf->Add("kr0.root");
  etf->Add("ch.root");
  etf->Add("c.root");
  TFile *f = new TFile("ch.root");
  cout << "Totally " << etf->GetEntriesFast() << " events" << endl;
  TCut cut = "1 "
    "&& x_f1>=-19.&&x_f1<=-11. "
    "&& pdcTaUWire[1][0]>=-15.&&pdcTaUWire[1][0]<=-1. "
    // "&& pid.betaz[0]()>0.6&&pid.betaz[0]()<0.7 "
    // "&& fabs(globalTracker.dx2())<4. " // good pid
    // "&& pid.poz()>1300. "
    "&& sqrt((pTAU[0]+7.)*(pTAU[0]+7.)+(pTAU[1]+1.5)*(pTAU[1]+1.5))<25. " // within TA
    "&& fabs(pT1[0]+6.2)<24. && fabs(pT1[1]+3.2)<24. " // within tstop
    "&& z0c>=19.5&&z0c<=20.5 "
    "&& (ca21 || ca22 || ca23 || ca24) "
    // "&& tof1>=135.&&tof1<=155. "
  ;
  const int N = 4;
  const char *az[N] = {"ca21", "ca22", "ca23", "ca24"};
  cout << "\033[36;1m" << cut << "\033[0m" << endl;

  etf->Draw(">>el", cut, "entrylist");
  TEntryList *el = (TEntryList *)gDirectory->Get("el");

  etf->Draw("z0c:aoz0c>>h1(1000, 1.9, 2.4, 1000, 14., 39.)");
  etf->SetEntryList(el);
  etf->Draw("et.tof2:aoz0c>>h2(100, 1.9, 2.4, 150, 60., 80.)");
  // etf->Draw("pid.aoz1():aoz0c>>h3(100, 1.9, 2.4, 120, 1.9, 2.6)");
  etf->Draw("pid.aoz1()*20.>>h3(75, 39, 49)");
  // etf->Draw("et.tof2:et.tof1>>h2(200, 138., 148., 200, 61., 74.)");
  // etf->Draw("pid.beta():pid.betaz[0]()>>h3(200, 0.59, 0.62, 200, 0.51, 0.64)");
  // etf->Draw("postMagTrackerX.GetHitPos(\"tofwB\"):aoz0c>>h4(200, 1.9, 2.4, 200, -700., 700.)");
  // etf->Draw("pid.aoz1():aoz0c>>h4(100, 1.9, 2.4, 120, 1.9, 2.6)");
  etf->Draw("pid.bg0():pid.bg()>>h4(500, 0.5, 0.9, 500, 0.5, 0.9)", "aoz0c>0.", "col");
  // etf->Draw("postMagTrackerX.ang()-globalTracker.GetPostTaTrack().ang():aoz0c>>h4(200, 1.9, 2.4, 200, 10., 20.)");
  // etf->Draw("postMagTrackerX.GetHitPos(\"tofwB\")>>h6(300, -700., 700.)");
  etf->Draw("pid.poz():pid.beta()>>h5(400, 0.5, 0.9, 400, 1100., 1600.)");
  etf->Draw("pid.aoz1()*20.>>h6(75, 39, 49)");

  TCanvas *c = new TCanvas("c", "ETF PID Analysis", 1680, 1050);
  c->Divide(3,2); char h[64]; TH1 *b;
  for(int i = 0; i < 6; i++){
    sprintf(h, "h%d", i+1); f->GetObject(h,b);
    c->cd(i+1)->SetGrid();
    if(2 != i) b->Draw("col");
    else b->Draw("e");
  } // end for

  c->cd(1);
  etf->Draw("z0c:aoz0c>>h1_0(1000, 1.9, 2.4, 1000, 14., 39.)", "", "same");

  char cmd[128]{};
  c->cd(6);
  for(int i = 0; i < N; i++){
    sprintf(cmd, "pid.aoz1()*20.>>ha%d(75, 39, 49)", i);
    etf->Draw(cmd, az[i], "same");
  } // end for over i
  double aoz[5]{};
  double par[12]{}; // gaussian for the partial aoz-s
  for(int i = 0; i < N; i++){
    sprintf(h, "ha%d", i); f->GetObject(h, b);
    b->SetFillColor(i+1);
    b->Fit("gaus");
    TF1 *gaus = b->GetFunction("gaus");
    par[3*i] = gaus->GetParameter(0); // constant
    aoz[i] = par[3*i+1] = gaus->GetParameter(1); // mean
    par[3*i + 2] = gaus->GetParameter(2); // sigma
  } // end for over i
  // etf->Draw("pid.aoz1()>>ha5(103, 1.9, 2.6)", az[5], "same");
  // ((TH1F *)f->Get("ha5"))->SetFillColor(5);
  c->cd(3);
  TF1 *aoz_f1[N];
  int color = 1;
  for(int i = 0; i < N; i++){
    sprintf(h, "aoz_f1%d", i);
    aoz_f1[i] = new TF1(h, "gaus", 39, 49);
    aoz_f1[i]->SetNpx(3000);
    color++;
    while(1 == color || 4 == color || 5 == color) color++;
    // aoz_f1[i]->SetLineStyle(2);
    aoz_f1[i]->SetLineColor(color);
  } // end for over i
  TF1 *aoz_f1_merge = new TF1("aoz_f1_merge", "gaus(0)+gaus(3)+gaus(6)+gaus(9)", 39, 49);
  aoz_f1_merge->SetNpx(6000);
  aoz_f1_merge->SetLineColor(4); aoz_f1_merge->SetLineWidth(4);
  aoz_f1_merge->SetParameters(par);
  // for(int i = N*3; i--;) aoz_f1_merge->FixParameter(i, par[i]);
  // for(int i = N; i--;) aoz_f1_merge->FixParameter(3*i+2, par[3*i+2]);
  for(int i = N; i--;) aoz_f1_merge->SetParLimits(3*i+2, 0.2, 0.4);
  ((TH1F *)f->Get("h3"))->Fit(aoz_f1_merge);
  aoz_f1_merge->GetParameters(par);
  for(int i = 0; i < N; i++){
    aoz_f1[i]->SetParameters(&par[3*i]);
    aoz_f1[i]->Draw("same");
  } // end for over i


  c->cd(5);
  for(int i = 0; i < N; i++){
    sprintf(cmd, "pid.poz():pid.bg()>>hp%d(400, 0.5, 0.9, 400, 1100., 1600.)", i);
    etf->Draw(cmd, az[i]);
  } // end for over i
  for(int i = 0; i < N; i++){
    sprintf(h, "hp%d", i); f->GetObject(h, b);
    b->SetMarkerStyle(3); b->SetMarkerColor(i+1);
    b->Draw("same");
  } // end for
  TF1 *faz[5];
  const double zz = 18.;
  for(int i = 0; i < N; i++){
    sprintf(h, "pv%d", i);
    // sprintf(h, "pv%d", i);
    // aoz[i] = (zz*2+1+i) / zz;
    faz[i] = new TF1(h, pbg, 0.4, 0.9, 1, 1);
    faz[i]->SetParameter(0, (aoz[i]/20.-0.494146)/0.738185);
    faz[i]->SetLineColor(i+1);
    faz[i]->Draw("same");
  } // end for

  c->Print("pt20m.pdf");

  TFile *fw = new TFile("pid20.root", "recreate");
  for(int i = 0; i < N; i++){
    aoz_f1[i]->Write("", TObject::kOverwrite);
  } // end for over i
  aoz_f1_merge->Write("", TObject::kOverwrite);
  ((TH1F *)f->Get("h3"))->Write("", TObject::kOverwrite);
  c->cd(3)->Write("", TObject::kOverwrite);
  fw->Close();
} // end pt

// p[0]: aoz; x[0]: beta
// return p/Z = poz(beta, aoz);
double pv(double *x, double *p){
  static const double cc = ETFMath::CC(), c0 = ETFMath::c0();
  const double aoz = p[0], beta = x[0];
  const double brho = aoz/cc*ETFMath::BetaGamma(beta);
  const double poz = brho * c0;
  return poz;
} // end function pv
double pbg(double *x, double *p){
  static const double cc = ETFMath::CC(), c0 = ETFMath::c0();
  const double aoz = p[0], bg = x[0];
  const double brho = aoz/cc*bg;
  const double poz = brho * c0;
  return poz;
} // end function pv
